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SUMMARY 

We present a high-order shifted Gegenbauer pseudospectral method (SGPM) to solve numerically the second-order one¬ 
dimensional hyperbolic telegraph equation provided with some initial and Dirichlet boundary conditions. The framework of the 
numerical scheme involves the recast of the problem into its integral formulation followed by its discretization into a system 
of well-conditioned linear algebraic equations. The integral operators are numerically approximated using some novel shifted 
Gegenbauer operational matrices of integration. We derive the error formula of the associated numerical quadratures. We also 
present a method to optimize the constructed operational matrix of integration by minimizing the associated quadrature error 
in some optimality sense. We study the error bounds and convergence of the optimal shifted Gegenbauer operational matrix 
of integration. Moreover, we construct the relation between the operational matrices of integration of the shifted Gegenbauer 
polynomials and standard Gegenbauer polynomials. We derive the global collocation matrix of the SGPM, and construct 
an efficient computational algorithm for the solution of the collocation equations. We present a study on the computational 
cost of the developed computational algorithm, and a rigorous convergence and error analysis of the introduced method. 
Four numerical test examples have been carried out in order to verify the effectiveness, the accuracy, and the exponential 
convergence of the method. The SGPM is a robust technique, which can be extended to solve a wide range of problems arising 
in numerous applications. 
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1. INTRODUCTION 

Second-order hyperbolic partial differential equations (PDEs) have been studied for many decades, as they 
frequently arise in many applications like seismology, acoustics, general relativity, oceanography, electromagnetics, 
electrodynamics, thermoelasticity, thermodynamics of thermal waves, fluid dynamics, reaction-diffusion processes, 
materials science, geophysics, biological systems, ecology, and a host of other important areas; cf. [1, 2, 3, 4, 5]. The 
range and significance of their applications manifest the demand for achieving higher-order numerical approximations 
using robust and efficient numerical schemes. In the present work, we establish a high-order numerical approximation 
to the solution of the following second-order one-dimensional hyperbolic telegraph equation: 

Utt + PiUt-f I52U = Uxx + 0<X<1, 0<t<T, (1-la) 


provided with the initial conditions given by 

(1.1b) 
(1.1c) 


u{x, 0) = gi{x), 0 < X < I, 
ut{x, 0) = g 2 {x), 0 < X < I, 
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and the following Dirichlet boundary conditions 

u{0,t) = hi{t), 0 < t < T, (1-ld) 

u{l,t) = h 2 {t), 0 < f < r, (1-le) 

where u is the unknown solution function, / is a given integrable function; gi,g 2 , hi, and /12 are some given functions; 
Pi and P 2 are some known constant coefficients. We shall refer to Equation (1.1a) provided with Conditions (1.1b) 
- (l.le) by Problem V. In fact, the telegraph equation (1.1a) models an infinitesimal piece of a telegraph wire as an 
electrical circuit, and it describes the voltage and current in a double conductor with distance x and time t [6]. The 
telegraph equation is in particular important as it is commonly used in the study and modeling of signal analysis 
for transmission and propagation of electrical signals in a cable transmission line [7, 8], and in reaction diffusion 
occurring in many branches of sciences [9, 10]. 

The numerical solution of second order hyperbolic PDEs has been studied extensively by a variety of techniques 
such as the finite element methods [11, 12], finite-difference schemes [13, 14, 15, 3], combined finite difference 

scheme and Haar wavelets [6], discrete eigenfunctions method [7], Eegendre multiwavelet approximations [16], the 

singular dynamic method [17], interpolating scaling functions [18], cubic and quartic B-spline collocation methods 
[19, 20], non-polynomial spline methods [21], the reduced differential transform method [22], etc. In the present 
work, we present a shifted Gegenbauer pseudospectral method (SGPM) for the solution of Problem V. The numerical 
scheme exploits the stability and the well-conditioning of the numerical integral operators, and collocates the integral 
formulation of Problem V in the physical (nodal) space using some novel operational matrices of integration 
(also called integration matrices) based on shifted Gegenbauer polynomials. The proposed method leads to well- 
conditioned linear system of algebraic equations, which can be solved efficiently using standard linear system solvers. 
The rapid convergence, economy in calculations, memory minimization, and the simplicity in programming and 
application are some of the features enjoyed by the present method. The current work is an extension to the works 
of Elgindy and Smith-Miles [23] and Elgindy and Smith-Miles [24] to second-order hyperbolic PDEs using shifted 
Gegenbauer polynomials. 

The rest of the article is organized as follows: In Section 2, we give some basic preliminaries relevant to Gegenbauer 
and shifted Gegenbauer polynomials. In Section 3, we derive the Eagrange form of the shifted Gegenbauer 
interpolation at the shifted Gegenbauer-Gauss (SGG) nodes. In Section 4, we derive the shifted Gegenbauer 
integration matrix and its associated quadrature error formula in Section 4.1. In Section 4.2, we construct an optimal 
shifted Gegenbauer integration matrix in some optimality measure, and analyze its associated quadrature error in 
Section 4.3. Section 4.4 gives the error bounds of the optimal shifted Gegenbauer quadrature. Section 4.5 presents 
the relation between the integration matrices of the shifted Gegenbauer polynomials and standard Gegenbauer 
polynomials. Section 5 introduces the SGPM for the efficient numerical solution of Problem V. Section 5.1 is devoted 
for the constructions of the global collocation matrix and the right hand side of the collocation equations. Section 5.2 
establishes the global approximate interpolant over the whole solution domain. Section 5.3 is devoted for the study 
of the convergence and error analysis of the proposed method. Pour numerical test examples are studied in Section 
6 to assess the efficiency and accuracy of the numerical scheme. We provide some concluding remarks in Section 8 
followed by possible future work in Section 7. Pinally, Appendix A establishes an efficient computational algorithm 
for the constructions of the global collocation matrix and the right hand side of the collocation equations. 


2. PREEIMINARIES 


In this section, we present some preliminary properties of the Gegenbauer polynomials and the shifted Gegenbauer 
polynomials defined on one and two dimensions. Moreover, we present the discrete inner product of any two functions 
for the shifted Gegenbauer approximations. 

The Gegenbauer polynomial c!^\x), of degree n G Z+, and associated with the parameter a > —1/2, is a real¬ 
valued function, which appears as an eigensolution to a singular Sturm-Eiouville problem in the finite domain [—1,1] 
[25]. It is a Jacobi polynomial, p!pp‘'^\ with a = P, and can be standardized so that: 




n!r(a-h 1) 
r(n + a-\- 




n = 0,l,2,.... 


( 2 . 1 ) 


Therefore, we recover the nth-degree Chebyshev polynomial of the first kind, Tn{x), and the nth-degree Eegendre 
polynomial, Ln{x), for a = 0 and 0.5, respectively. The Gegenbauer polynomials can be generated by the following 
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three-term recurrence equation: 

(n +2a)C^‘^^(x) = 2(n + a)xC^°‘^x) - n = 1,2,3,..., (2.2a) 

starting with the following two equations: 

C(“^(x) = l, (2.2b) 

C['"\x) = x, (2.2c) 

or in terms of the hypergeometric functions, 

= 2 F 1 (-n,2a + n;a+ , n = 0,1, 2,..., (2.3) 


where 2 F 1 {a,b-,c-,x ), is the first hypergeometric function (Gauss’s hypergeometric function), which converges if 
c ^ Z“ U {0}, for all |a:| < 1, or at the endpoints x = ±1, if Re[c — a — b] >0. The leading coefficients of the 
Gegenbauer polynomials C^\x), are denoted by k!^\ and are given by the following relation: 


^(a) ^ nn -1 + o) r(2o: + 1) 

” r(n + 2a) r(a + l)’ 


71 = 0 , 1 , 2 ,.... 


(2.4) 


The weight function for the Gegenbauer polynomials is the even function (a:) = (1 — The Gegenbauer 

polynomials form a complete orthogonal basis polynomials in [—1,1], and their orthogonality relation is given 
by the following weighted inner product: 




Cl^\x)Cl^\x)w^^\x)dx 



— Gr7 


(2.5) 


where 




2i-2“7rr(n + 2a) 
n\(n +a)T"^{a) 


( 2 . 6 ) 


is the normalization factor, and 5m,n is the Kronecker delta function. We denote the zeroes of the Gegenbauer 
polynomial {x) (also called Gegenbauer-Gauss nodes) by and denote their set by . We 

also denote their corresponding Chris toff el numbers by 07^“^, /c = 0,..., n, and define them by the following relation: 


(^S) ^ = 0,l,2,...,n. 


(2.7) 


3=0 


Throughout the paper, we shall refer to the Gegenbauer polynomials by those constrained by standardization (2.1). 

Let L be some positive real number. The shifted Gegenbauer polynomial of degree n on the interval [0,L], 
is defined by = clt\2.x/L— 1). The shifted Gegenbauer polynomials form a complete 1,1]- 

’ '^L 

orthogonal system with respect to the weight function. 






uu^ {x) = {Lx — X^Y 

and their orthogonality relation is defined by the following weighted inner product: 


Jo 


c 


(«) 


A — A 

(a) ^FTl^n — 


( 2 . 8 ) 


(2.9) 


where 


A 


(a) ^ 
L,n 



( 2 . 10 ) 


is the normalization factor. For a = 0 and 0.5, we recover the shifted Chebyshev polynomials of the first kind and the 
shifted Legendre polynomials, respectively. We denote the zeroes of the shifted Gegenbauer polynomial 
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(SGG nodes) by /c = 0,..., n, and denote their set by We also denote their corresponding Christoffel 

numbers by f,,k = Q,... ,n. Clearly 


“i.fc = \ (4“1 + l) . fc = 0,..., n; 


w 


(a) 


2a 


= (n) k = 0,...,n. 


''L^n,k \ 2 ' 

If we denote by P„, the space of all polynomials of degree at most n, n € Z+, then for any </> e P 2 n+i. 


( 2 . 11 ) 

( 2 . 12 ) 


^ ^(x)4">Mdx= Q + „.<“)M<ix= (§) 

n 


2 V "4 


0 ) 


(2.13) 


using the standard Gegenbauer-Gauss quadrature. With the quadrature rule, we can define the discrete inner product 
(■; Oi.n’ of two functions u{x) and v{x) defined on [0, L], for the shifted Gegenbauer approximations as follows: 


n 

(m, v)L,n = ^ ^ j) ^ ■ 


(2.14) 


j=0 


In two dimensions, we can define the bivariate shifted Gegenbauer polynomials by the following definition: 
Definition 2.1 

Let Ici^hx)} , be a sequence of shifted Gegenbauer polynomials on Dl = [0,L]. The bivariate shifted 

L ’ J n=0 ^ 

Gegenbauer polynomials, | i.rC^m^x, f) I , are defined as 

L ) n,m=0 

i.rC^^l{x, t) = (x, t) G Dl, = [0, 1] X [0, r], f, r G M+. (2.15) 

The family , forms a complete basis for L'^{Df.^). They are orthogonal on with 

respect to the weight function (x, t) = (x) (t), since 

[ [ i,TCl^}n{x,t)i,TC^^^{x,t)w\‘^}{x,t)dxdt= f Cl“j{x)cl‘l\x)wl‘^\x)dx- f dt 

Jo Jo Jo Jo 

(2.16) 

(2.17) 

(2.18) 


{n,m) = {s,k), 

\ 0, O.W., 


where 


, Vr? m 


In the next section, we highlight the modal and nodal orthogonal shifted Gegenbauer interpolation, and derive the 
Lagrange form of the shifted Gegenbauer interpolation at the SGG nodes. 


3. ORTHOGONAL SHIFTED GEGENBAUER INTERPOEATION 


The function 


Pnfix)=Y.fjC'tj{x), 
j=0 


-((a)/ 


(3.1) 


This is a preprint copy that has been accepted for publication in Numerical Methods for Partial Differential Equations. 



HIGH-ORDER NUMERICAL SOLUTION OE THE TELEGRAPH EQUATION 


5 


is the shifted Gegenbauer interpolant of a real function / defined on [0, L], if we compute the coefficients fj so that 


Pnf{xk) = f{xk), k = 0,...,n, 


(3.2) 


for some nodes Xk G [0,L],k = 0,... ,n. If we choose the interpolation points Xk to be the SGG nodes, then we 
can simply compute the discrete coefficients using the discrete inner product created from the Gegenbauer-Gauss 
quadrature by the following formula: 


fj = 




L,j>L,n 


a 


(a) 

L,J 




c 


(a) 

LJ 


X 






(a) 

L,j k=0 


Aa 
.,k JL'. 


r'(“) 
,fc J 


a) \ 

'.,n,k J ’ 


J = 0, 


(3.3) 


where k ~ f Equation (3.3) gives the discrete shifted Gegenbauer transform. To construct the 

shifted Gegenbauer integration matrix, we need to represent the orthogonal shifted Gegenbauer approximation as 
an interpolant through a set of node values (nodal approximation) instead of the modal approximation given by 
Equation (3.1). Substituting Equation (3.3) into (3.1) yields 


Pnf{x) 


^(a) 

j=0 ^L,j k=0 


(a) Aa) ^{a) (^(. 
^L,n,k J L,n,k ^L,j 




E 

fc =0 . 


n 




lL,n,k- 


(3.4) 


Hence the Eagrange form of the shifted Gegenbauer interpolation of / at the SGG nodes can be written as: 

n 

Pnfix) = 

fc=0 


(3.5) 


where „ k(^)’ Eagrange interpolating polynomials defined by 


lb 

(^S) C'S (2^). k = 0,...,n. 

j=o 


(3.6) 


Theorem 3.1 

The functions ^(a;), fc = 0,..., n, defined by Equation (3.6) are the Eagrange interpolating polynomials of the 
real-valued function / constructed through shifted Gegenbauer interpolation at the SGG nodes. 


Proof 

To show that ^(t), A: = 0,..., n, are indeed the Eagrange interpolating polynomials of the real-valued function 
/, we need only to show that ki^^Ln t) = i,k = 0,... ,n. Since 


/•(“) 

^L,n,k 


{x) 


2a 

2 I tn. 




lb .■ 

SE (")“’)■ Q“'(43)crM 


j=o 




SE("r) 


j=0 


r (j -I- a -I- 2) 




S) 


zu 


lb .< 

cr (xS) cf>M, 


J=0 


(3.7) 
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where (7(a;), are the Gegenbauer polynomials standardized by Szego [25], and 


A<“> = 


c: 


(a) 


2 ^ 22 “-ir 2 (j + a+i) 


!,(<») (j + a)j\T{j + 2a) 

By Christoffel-Darboux Theorem (see [26, Theorem 4.4]), 

-(2a+i) 2) r2(2a + 1) ^^n+i{x) 




(a) 2 


r(n + 2 a + l)r2(a + l) 


X — X, 


(a) 


(3.8) 


(3.9) 


Hence = 0 Vi 7 ^ fc, since ) = 0, i = 0,..., n. For i = k, and using L’HopitaTs rule, we find 

that 


i(“) /.„(“) 


C 


‘) A _ _(“) 2 r(n + 2) r^(2a; + 1) - (a) x a(q) / (a) n 

.n.k\-^n.k) ~ ‘^n.k , o„, , iAt^2A„, , 1 A y'^n,k> ^^'-'n+l\-^n,k)- 


-L,n,kK n,kJ n,k r{n + 2a + l)T^ia + 1) 

Using Equation (5.19) in [26], we can easily show that 


= (n + 2 a + 1 ) (^(^1 - ( 4 “l) ) 


=1. k = 0,...,n. (3.11) 


Hence 




(3.10) 


£' 


/t _-1 

«) >(“)("..(“) ) 

'.,n,ky L,n,i) ~ ^L,n,k \ J \^L,n,kJ ^L,n,i) 


3=0 


' I' _ 1 

”2E (V’)' V"’ (VI) c<"’(V“>) = *,*■ 


3=0 


(3.12) 

□ 


4. THE SHIFTED GEGENBAUER INTEGRATION MATRIX 

Suppose that a real-valued function / is approximated by the shifted Gegenbauer interpolant P„/ given by Equation 
(3.5). The shifted Gegenbauer integration matrix calculated at the SGG nodes is simply a linear map, which 
takes a vector of {n + 1) function values, F = (/l,oi /l, 1 i ■ ■ ■ j fL,n)'^, to a vector of (n -|-1) integral values 


ll“) = ( / Pnf{x) dx, / Pnf{x) dx,..., / Pnf(x) dx j , 
Jo Jo Jo / 


such that 


l(“) = F. 


(4.1) 


The integration matrix = {p^l\ fc), 0 < i, fc < n, is the first-order square shifted Gegenbauer integration matrix 

of size (n -I- 1), and its elements, p^j^\ = p^j^\ ^.(a), 0 < i, A: < n, can be constructed by integrating Equation (3.5) 

on [ 0 , L], such that 


/ 


Pnf(x)dx = J^fi‘^l^ / 6^lj^{x)dx 

k=0 


E 41.E(^S)' ^S(4E) 

fc =0 V J =0 ^ ^ -^0 




_y-Ai) Aa) 

~ ^^L,i,k J L,n,k' 


(4.2) 
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Hence 


»( 1 ) (q 

PL,^,k = 


Ifc (^S) / C'ij (^) dx, i, fc = 0,..., n. (4.3) 

i=o 


We refer to the shifted Gegenbauer integration matrix, and its associated quadrature (4.2) by the S-matrix and 
the S-quadrature, respectively. 


4.1. Error Analysis of the S-Quadrature 

The following theorem highlights the truncation error of the shifted Gegenbauer quadrature associated with the shifted 
Gegenbauer integration matrix 

Theorem 4.1 


Let /(x) G [0, L], be interpolated by the shifted Gegenbauer polynomials at the SGG nodes, ^ G then 


j(«) 


there exist a matrix P^^^ = 0 < i,j < n, and some numbers J G (0, L), i = 0,... ,n, satisfying 


(q) 


L,i,j 


L,n,iy 


f ’ ’ f{x)dx = 

•^0 fc =0 


(4.4) 


where pp \ i, fc = 0,..., n, are the elements of the matrix P^^\ defined by Equation (4.3), and 


(x(“) 


L.n.i ’ 




2 / {n + l)\Ki%Jo 


f 




C'L,n-ri(^) dx. 


Proof 

Set the error term of the shifted Gegenbauer interpolation as 


Rn[x) = /(x) - Pnf{x), 


and construct the auxiliary function 


Y{t) = - 


Rn{x) 


ctUM 


ctUY)- 


Since / G L], and Pnf G L], it follows that Y G G”+^[0, L]. For t = xY p have 


= RnYpY)- 


Rn{x) ^(g) 


Yl.L) = 0 , 


'^L,n,iy ^v‘^L,n,z/ ( \ L.,rL,i 

'L.n+iy^) 


since x^“^ are zeroes of i?„(x). Moreover, 


r{x) = i?„(x) - = 0. 

Ch,UY) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 


(4.9) 


Thus Y G [0, L], and H is zero at the (n + 2) distinct nodes x, x^^ i = 0,..., n. By the generalized Rolle’s 
Theorem, there exists a number ^ in (0, L) such that y (^) = 0. Therefore, 


0 _ y (n+l) (-f \ _ D(n+1) ("f \__ d ^{a) , ,x 




L,n+1 

n+1 










" " U; ^ ^ ”^^gk^,(x)- 




(4.10) 
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Since P„/ e P„, (x) is identically zero, and we have 


0 = ) (n + iy.K, 

^+1 y*(n+l)('^^ 


(a) -Sn(^) 


■■■ f(x) = Pnf(x) + ( ^ 


(n + iy.K7, 






Jo) 


r-fixux^r-i^nxux. r 

io Jo V2/ (n + l)!it:i“\ Jo 

n 

= Vj5^^^ .(«) 

/ u ^L.i.kJ L,n,k ^ ^L,n y^L,n,i'’ y * 




fe =0 


(4.11) 


(4.12) 


(4.13) 

(4.14) 

□ 


4.2. Optimal S-Quadrature 

To construct an optimal S-quadrature to approximate the definite integration f{x) dx, of an integrable function 
/, for any arbitrary integration node Xi S [0, L], we follow the idea presented by Elgindy and Smith-Miles [23], and 
seek to determine the optimal Gegenbauer parameter a*, which minimizes the magnitude of the quadrature error 
E^yy\{xi,^), at each integration node Xi,i = 0, ■ ■ ■ ,n. Define the smooth function 77L,i,n(a), such that 


VL,i,n{a) 


K, 


n+1 


f 




dx. 


(4.15) 


The values of the optimal Gegenbauer parameters a*, can be determined through the following one-dimensional 
optimization problems: 

Find a* = argmin 77 ! , „(a), i = 0,. .. ,n. (4.16) 

a>-l/2 

Problems (4.16) can be further converted into unconstrained one-dimensional minimization problems using the 
following change of variable: 

Oi = (l-t-2G;) 1. (4.17) 

We refer to the optimal shifted Gegenbauer integration matrix and its associated quadrature established through the 
solution of Problems (4.16) by the optimal S-matrix and the optimal S-quadrature, respectively. Notice here that for 
each integration node Xi, an optimal Gegenbauer parameter a* is determined, and the optimal S-quadrature seeks 
a new set of SGG nodes as the optimal shifted Gegenbauer interpolation nodes set corresponding to the integration 
node Xi. We denote these optimal SGG interpolation nodes by k = 0,... ,m, for some m € Z+, and we 

call them the adjoint SGG nodes, since their role is similar to the role of the adjoint Gegenbauer-Gauss nodes 

(a*) 

’ j,, fc = 0, ..., m, constructed in [23]. Notice also that the choice of the positive integer number m is free, which 
renders the optimal S-matrix a rectangular matrix of size (n -I-1) x (m -I-1) rather than a square matrix of size 
(n -I- 1), as is typically the case with the standard S-matrix. Denote the optimal S-matrix by = {p^l\ f.)0 = 
0,..., n; fc = 0,..., TO, where p^j^\ = p^l\ are the matrix elements of the ith row obtained using the optimal 

value a*. The definite integral f{x) dx, is then approximated by the optimal S-quadrature as follows: 



m 

dx « 


/c =0 




i = 0 ,... ,n, 


(4.18) 


where 


,k — J \^L,mA 


A, i = 0,...,n;k = 0, 


, m. 
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4.3. Error Analysis of the Optimal S-Quadrature 

The following theorem describes the construction of the optimal S-matrix elements, and highlights the truncation 
error of the associated optimal S-quadrature. 

Theorem 4.2 
Let 


SL,u,m = = 0, * = 0,..., n; fc = 0,..., m}, L e M+, 


n, m 


e z+, 


(4.19) 


be the adjoint set of SGG nodes, where a* are the optimal Gegenbauer parameters in the sense that 


a* = argmin77|_^_^(a). 

q>-1/2 


(4.20) 


Moreover, let f{x) G be a real-valued function approximated by the shifted Gegenbauer polynomials 

expansion series such that the shifted Gegenbauer coefficients are computed by interpolating the function f{x) at 

(a*) 

the adjoint SGG nodes mik ^ SL,n,rm i = 0,. ■. ,n;k = 0,... ,m. Then for any arbitrary integration nodes Xi G 

[0, L],i = 0,... ,n, there exist a matrix = {p^l\ j),i = 0,... ,n;j = 0,... ,m, and some numbers G [0, L], 
satisfying 


pXi ''t ^ 

/ f{x)dx = Y^ ft.il.i.k + - ^ 0 , 

k=0 


where 


PL,i,k — ^L,m,k 2-^ \^^L,j J ^L,j \^L,m,i,k) / ^L,j 
j=0 "'o 




"^+1 /(™+l)(^^) 


2™ (m -L 1)! 



(4.21) 

Epj \x) dx; 

(4.22) 


(4.23) 


Proof 

The function 


Pmfix)=Y.kjC^L] 

3=0 


(4.24) 


is the shifted Gegenbauer interpolant of the real function / defined on [0, L], if we compute the coefficients fj so 
that 

(4.25) 


Pmf{z^L,ii,k) = f = 0,..., n; fc = 0,..., m. 


Hence the discrete shifted Gegenbauer transform is 

^ L.m f 


fi,k — 




c 


(«•) 

L,3 




c 


L,3 


A 


1) 


(a) 


L.j k—{) 


AO \ ,■ 

Z_^^L,m,kJ L,m,i,k'~’L,j j > ‘ 


= 0,...,n; k = 0,... ,m. 


(4.26) 

Following the approach presented in Section 3, we can easily show that the Lagrange form of the shifted Gegenbauer 
interpolation of / at the adjoint SGG nodes can be written as: 


Pmf{x) — fL.fi.i.k PL.ln.i.kjx), 


(4.27) 


fc=0 


where Cpp ^ ki^), are the Lagrange interpolating polynomials defined by 


m _ ^ 

Gi^,k(x) = {a^,k) Cppix), i = 0,...,n;k = 0,...,m. (4.28) 


3=0 
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Therefore, 



Pmf{x) dx 





m 

/ j L,m,i^k I 

1 _n 0 








TTl 1 nQ[^ 

V«) ^ / 

^L,m,k 2_^ y^L,j J ^L,j y^L,m,i,kJ 

3=0 -^0 


ci“j\x)dx 


f(“*) 

L,m,i,k 


y- (1) .K*) 

/ V ^L,i^kd L,m,i,k' 


fc=0 


^PL,i,kJL,m,i,k ^ ^L,m 


fc=0 


where the quadrature error term, (x^, ^i), follows directly from Theorem 4.1 on substituting the value of a with 
a*, and expanding the shifted Gegenbauer expansion series up to the (m + l)th-term. □ 


4.4. Error Bounds of the Optimal S-Quadrature 

To smdy the error bounds of the optimal S-quadrature, we require the following two lemmas. 
Lemma 4.1 

The maximum value of the shifted Gegenbauer polynomials G^“^(a:), is given by 



ctA) 




n\T{2a) 

f f+«-l^ 

max 

x^[0, L] 

— 


= < 

L“[0. L] 

r(n + 2a) 

1 t ) 


1, n > 0 A a > 0, 


n „ I 1 
- G Z+ A -- < a < 0, 




n^^lA—- <a<0, 


(4.29a) 

(4.29b) 

(4.29c) 


where Zg = Z+ U {0}, is the set of all non-negative integers, and > 1, is a constant dependent on a, but 
independent of n. Moreover, for odd n > 0, and — ^ < a < 0, the maximum value of C^^\{x), is bounded by the 
following inequality 


max 

XG [0, L] 




C 


(a) 

L,n 


< 

L”[0.L] 


2 n!r( 2 a) |a| 

yjn{2a + n) r(n + 2a) 



(4.30) 


Proof 

The proof is straightforward. Indeed, Equalities (4.29) follow using Equation (2.1), Eemma 2.1 in [24], and Equation 
(7.33.2) in [25]. Since for a certain value of 0 < x < L, monotonically decreases for increasing values of 

a in the range —1/2 < a < 0, as n —?> oo; cf. [23, Appendix D], then 




C 


(a) 


L“[0,L] 


is also monotonically decreasing for increasing values of a G (—1/2,0). Since 


lim = 1, 

a—>-0 


(4.31) 


(4.32) 


then A^°‘l > IVa G (—1/2,0). Einally, Inequality (4.30) follows using Equation (2.1) and Equation (7.33.3) in 
[25]. □ 


Lemma 4.2 

Eor a fixed a > —1/2, the factor (n -I- l)!iT^“\, is of order n^/^““(2n/e)”, for large values of n. 
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Proof 

The lemma is a more accurate version of Lemma 2.2 in [24] by realizing that (n + 1)! = (n + 1) • n! « n • 
a/2 TT n [nje)^, as n —>■ oo. □ 


The following theorem gives the error bounds of the optimal S-quadrature. 


Theorem 4.3 (Error bounds) 

Assume that /(x) S C""+^[0,L], and < A G M+, for some number to G Zq . Moreover, let 

Jo' f{^) dx, be approximated by the optimal S-quadrature (4.18) up to the (to + l)th shifted Gegenbauer expansion 
term, for each integration node x^,i = 0,..., n. Then there exist some positive constants and D^'\ 

independent of to such that the truncation error of the optimal S-quadrature, is bounded by the following 

inequalities: 


E 


(a*) 

L,m 


^ 2 - 2 m-ip (q, XiL'^+^T {m + 2a + 1) 
r (2a -L 1) r (to -L 2) T (to -L a -L 1) 


1, TO > 0 A a > 0, 

(m-|-l)!r( 2 Q) / + a — 

r(m-|-2a-ri) I rri+l 



m+1 


2 


eJ'P 

A2-2™-ir(a) |a| x*L’"+i 

( f+« ^ 


+(to -1- 1) (2a -1- TO -1- l)r (to -1- a -1- 1) 

1 f ) 


y GZ+A-^ <a<0, 


GZ+A-i<a<0 

(4.33) 

(4.34) 


for all i = 0, 


E 


(“*) 

L,Tn 


< 


B 


(a-) e™L™+i; 


22m+l y^m+3/2—a* 


/ 7-m+l 

(a*) L/ 


., n, where i?|“‘ ^ = AD\“'\ and ’ = B]“'’D 


B. 


(“*) 


Xi 


22771+1 771+3/2 ’ 


jK*) 


A“*), 


a* > 0 A TO » 1, 

-^ < a* < 0 A TO >> 1, 
(a,*) 


(4.35a) 

(4.35b) 


Proof 

The proof is straightforward using Equation (4.23), and Eemmas 4.1 and 4.2. □ 


4.5. The Relation Between the S-matrix and the Gegenbauer Integration Matrix 

Eet = (p/fe), / ^ = 0,..., n, be the first-order square Gegenbauer integration matrix of size (n -I-1), constructed 
by Theorem 2.1 in [23], and consider the Eagrange form of the shifted Gegenbauer interpolation of a real-valued 
function / at the SGG nodes, Pnf, given by Equation 3.5. Since 


" "■ r+i i T r 

= ^PLikiPnfJLl^k = ’ ’ Pnf{x)dx =- 

k=0 k=0 40 J- 

L 


-1 
L 


Pnf ( -^{x+1) ) dx 


Pnf (^{X+ 1)') dx = Pnf (xj^l + l)) = § '^P^k (Pn/JL^, 


o / ^ ^ik L,n,k^ 


(4.36) 


k—0 


where = Pnfix^L^ f.)^ then it follows that 



AD 

PL,i,k 

L .(1) 

= -^Pik , 

o' 

II 

(4.37) 

or in matrix form. 







■p(t) _ 

^:p(D 

2 

(4.38) 

Moreover, using the change of variable 






ti = 

P, 

—(si -L 1), 

i = 0,..., n, 

(4.39) 


This is a preprint copy that has been accepted for publication in Numerical Methods for Partial Differential Equations. 




















12 


KAREEM T. ELGINDY 


and Cauchy’s formula for repeated integration, we can calculate the g-fold dehnite integrals of Pnf on [0, Jj * = 
0,..., n, as follows: 


(a) 


r^L,n,i 

= / / •■• / / Pnf{to)dtodti...dtq_2dtq-i 

1 -n J-n *^0 •/O •/0*/0 

L 


fc =0 fc =0 

^ p2tq_ilL-l i‘2t2!L-l i‘2ti!L-l 

2 


/-I 


'-1 


/-I 


/-I 


-Pn/ ( ^(so + 1) ) dsodsi . . . dSq-2dSq-l 


(f)’ P~^ ■ ■ •/" P Pnf pso + 1)) dsods, . ..dSq.2dSq.^ = P' /S, 


1 


2y (g-l)!7_ 


j,O’--./(f(-i) 


dt = — 


^ _ . k 

^ ^ Ai^O 

T \ Q 1 

iv \ 1 X—^ ^ 


2) (g-1)! f 




(1) (s<^) _ ^(“) 

^n,k 




i\OL} 

d L,n,k ’ 

(4.40) 


where 


(pl 1 , i, fc = 0,..., n, is the S-matrix of order q,q G Z+. Hence, 


s(9) 


S(9) 


, (“) ( 
]^\1 \ xl ' — 


PLs,k 1^2^ \^2y ( 9 -I)! 




9-1 


i,fc = 0, ...,n, 


(4.41) 


or in matrix form. 


p(9) = ( M p(9) = ( ^ 


2y (g- 


p^_ o Jl.n+l) - 


J„+i,i ) ) oP(l) 


(4.42) 


where = [x^no^^^ni^ ■ ■ ■ ja:i“n,]^, Ji,j is the all ones matrix of size i x j, “®” and “o” denote the Kronecker 
product and Hadamard product (entrywise product), respectively. Hence the S-matrices of higher-orders can be 
generated directly from the first-order Gegenbauer integration matrix. Similarly, we can show that the optimal S- 
matrices of distinct orders are related with the Gegenbauer integration matrices constructed by Theorem 2.2 in [23] 
by the following equations: 


„(?) 

PL,i,k 






T) \ 

,i,k j 


q-l 


(9-1)! 


Pi,k I 


i = 0,..., n; fc = 0,..., m. 


(4.43) 


The optimal S-matrices of distinct orders can therefore be calculated efficiently using Equations (4.43), and 
Algorithms 2.1 and 2.2 in [23]. Since the convergence properties of the S-quadratures and the optimal S-quadratures 
are inherited from the convergence properties of the corresponding Gegenbauer integration matrices and optimal 
Gegenbauer integration matrices, respectively, the following useful result is straightforward. 

Corollary 4.1 (Convergence of the optimal S-quadrature) 

Assume that f{x)GC"^^^[0,L], and < A e M+, for some number m S Z+. Moreover, let 

Jq" f{x) dx, be approximated by the optimal S-quadrature (4.18) up to the {m + l)th shifted Gegenbauer expansion 
term, for each integration node Xi G [0, L], f = 0,..., n. Then the optimal S-quadrature converges to the optimal 
shifted Chebyshev quadrature in the L°°-norm, as m —> 00 . 


Proof 

The corollary can be established easily using Equations (4.43) and Theorem 2.4 in [23]. □ 
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5. THESGPM 


We commence our numerical scheme by recasting Problem V into its integral formulation. Thus twice integrating 
Equation (1.1a) with respect to t yields, 


t n(72 


u{x^t) — K{x^t) + Pi / u{x,a)da + P2 

Jo 

where 

Using the substitution 


II 


t pG2 


u(x, tJi) daida2 = 


0 Jo 


K{x,t) = {Pit + l)gi{x) -\- tg2{x). 
Uxx{x,t) = pix^t), 


i^^XX (^5 f^i) + fix,cri)) daida2, 

(5.1) 

(5.2) 

(5.3) 


for some unknown function p, we can recover the unknown solution function u and its first-order partial derivative 
Ux in terms of p by successive integration as follows: 


Ux{x,t) = / p{a,t) da + Ci{t), 


u{x, t) = 


I 
III 


p{ai , t) dai da2 -\- xci {t) + C2 {t ), 


(5.4) 

(5.5) 


where ci{t) and C2{t) are some arbitrary functions in t. Using Dirichlet boundary conditions (l.ld) and (l.le), we 
find that 


ci(t) = j ( h2{t) - hi{t) - 


I pG2 


p{ai,t) dai da2 ; 


0 ^'O 


C 2 (t) = hi{t). 

Eet 9 i^x = x/l, and define the function p{x, t), such that 

p{x,t) = Oi^x {h2{t) - hi{t)) + hi{t), 

Moreover, let 


(5.6) 

(5.7) 

(5.8) 


nX pGq — l pO'2 


IqlIWx,t))= / .../ / p{ao,t)daodai ... 

Jo Jo Jo Jo 

d(Tq—2d(7q— i, 

(5.9) 

pt r<yq-l rG2 pCTl 

/ / ■•■ / / 0(x,fJo)<i67o<i(7i ... 

^0 Jo Jo Jo 

d<Jq_2d(7q_ , 

(5.10) 

and define the operator 

T r(^) p\ r(^) 

^l,x — J-2,x ' 


(5.11) 

Then we can simply write the unknown solution u{x, t) as follows: 



pX pG2 pi p(72 

u{x,t) = / / p{ai,t) daida 2 — di,x / / p{ai,t) dai da 2 + Pix,t) 

Jo Jo Jo Jo 

= Jl,x 

(5.12) 

Hence Equation (5.1) becomes 

ftp ='it, 


(5.13) 

where 



G = Jl,x + Pi l[*l Jl,. + P2 Jl,. - 4% 


(5.14) 

^ = P- (pil[^}+P2li*}) P + I^y-, 


(5.15) 
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Ipix, t) = k{x, t) — Ipix, t). 


(5.16) 


If we expand the unknown function (j){x,t), in a truncated series of bivariate shifted Gegenbauer polynomials as 
follows: 


Nt 




(5.17) 


n—0 m=0 


then we can compute the continuous coefficients of the truncation, n = 0,..., m = 0,..., Nt, using the 
two dimensional weighted inner product as follows: 


(j), l,TCn,m{x, t)^ jT 1 ,tC^%{x, t) (x, t) dx dt 

/o' fo (x, t) dx dt 




(5.18) 


I , T 


Instead, we lay a grid of SGG nodes, ^ = 0, ■ ■ ■, = 0,..., iVj, on the rectangular domain 

and approximate the function (j) by interpolation at those nodes. Let 


^s,k = <t> ,s = 0,...,N,-,k = 0,...,Nt. 


(5.19) 


The polynomial interpolant of cj) in two dimensions can be written in terms of the discrete coefficients (j>n,m, or in the 
equivalent Lagrange form as follows: 


Nt 


Wx Nt 


PN:„,Nt4'{x, t) — 'y ^ ^ ^ ^n,m l,TCn,m{X: t) — ^ ^ ^ ^ 4's,k 1,tP- ^)’ (5.20) 


n—0 m—0 


s—0 k—0 


where 


Clearly, 


^T,Nt,j'^ ~ i, s = 0,..., N^^j, fc = 0,..., Nt, 


(5.21) 

(5.22) 


and by construction, we find that 


Wx Nt 


[ PN,„,NMx,tl°^l^j)dx = '^'^<l)s,k [ 

"^0 s =0 fc =0 

= E^S, 


,(<») 

'<,Nx 


, (x ) dx 




i = 0,...,N;,;j = 0,...,Nt 


(5.23) 


s=0 

/ PNt,,Nt4>{x\2i^ t,t) dt = y/ y/ (j)s^k / l,rP^Nl,Nt,8,k^^\,,N,„,V^) 

^ 0 o — n 7 n 0 


S = 0 fc = 0 
Nt 

s(l) 


= '^PT.ik<l^k,i, Z = 0, . . . ,iVa;; j = 0, . . . ,7Vt. 


(5.24) 


fe =0 


To determine the bivariate discrete shifted Gegenbauer transform, we can first determine the intermediate values. 


Vn y‘'T,Nt,kJ ’ 


n = 0,..., N:c', k = 0,..., Nt, in the x-direction such that 


, ZVx 

“ 2 ^i,n , n = 0,..., iVa;; fc = 0,..., iVj. (5.25) 

^ ^ s=o 

kn ic) 
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Therefore, 


N, Nt 


l,T^n,m 


, 2 LT^Af?.Af,.s.fc^s,fc ^r.Nt.k) ^ U — 0 , . . . , N^', 171 — 0, - 


l,T 


s—0 k—0 


where 




(a) (a) 




Nu 

(5.26) 

(5.27) 


s = 0,... ,lVx; k = 0,...,Nt, 

are the two-dimensional Christoffel numbers corresponding to the SGG nodes, ■, , z = 0,..., Nx',j = 

0,..., Nt. To find the equations for the grid point values we require that (j) satisfies the integral formulation of 
the hyperbolic telegraph PDE (5.13) at the interior SGG nodes such that 


: i’ — Oj ■ ■ ■ ; j — 0; ■ ■ • ; 


(5.28) 


where 




-I' 


(t) 






At) 




i = 0,...,N,;j = 0,...,Nu (5.29) 

/ = 0,...,iVt, 

7 \ / 


V' 


+ /' 


(U 




lpi,j — Ip ’ * — 0, ... , iVr; j — 0, . . . , Nt- 


Let 


a;^“^ - Z- 


V’i.O.fc) = tp J.fc) > i = 0, ■ ■ ■, N^; j = 0, ■ ■ ■, Nt; k = 0,..., Mt, 


(5.30) 

(5.31) 

(5.32) 

(5.33) 


for some Mt € Z+. We can approximate the terms in Equation (5.28) using the S-quadrature and the optimal S- 
quadrature as follows: 


^ ~ ^ = ^Pi% - ^i,AcL. 


(2) rh 

iV:c+l,fc I 




fc=0 


r(0 


T,Nt,3 


{Ot) 


r(0 






Nt 


(5.34a) 

(5.34b) 


fc=0 

Nt 


■ N, 




(5.34c) 


fe =0 


^s=0 


s=0 


Nt 


■ N, 




3 .p^i^^t) = Y.Prik[Y.Pkls-^i,xn T.Pk 

K — O 


(2) ^ 


vs=0 


s=0 


r(^) 


Mt 

(c) G (i) « -ff \c) G t) = ^pi^k V'^,0■.fc)> 


-r,Nt,j 


-r,Nt,j 


r(0 


,(«) 


r(*) 






fc=0 

Mt 


tM 


tM 


t^t.j 


V' = '^PrikNiPk); 

k^O 
Mt 

f = Y.P?j,kk{3,k)- 


(5.34d) 

(5.34e) 

(5.34G 

(5.34g) 


fc=0 
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Hence the discrete analogue of Equations (5.28) can be written as: 


J, 




\ ^ ’ T,Nt,j \ ^ \ ’ T.Nj.j ’ ’ T.JVj.j ’ 


.(°) 






4>{x, t) = 

(5.35) 


or simply in shorthand notation as: 




= Lr-^ 


Mt 


j — l^T ^ 




(5.36) 


The solution of the linear system (5.36) provides the values of = 0,... ,Nx]j = 0,..., A^t, at the SGG nodes, 
G = 0,---,^x;j = 0,...,Nt; hence the discrete coefficients = 0, ■ ■ ■, N^;j = 0,...,Nt, 


from Equation (5.26). Notice that the coefficient matrix, of the linear system (5.36) generated by the 

SGPM is full, but in compensation, the high order of the basis functions gives high accuracy for given Nx,Nt;Mt. 

Using Equations (5.12) and (5.20), we can approximate the unknown solution u by the bivariate shifted Gegenbauer 
interpolant, as follows: 


Nt 


U{x,t) « PN^,NtU{x,t) = ^ X] i’n,mJl,x l,TC^%{x,t) +1p{x,t). 


(5.37) 


n=0 m—0 


Hence the approximate values of the unknown solution u can be determined at the SGG nodes, ^x\°‘^ j)’ 

i = 0,Nxij = 0,..., iVt, through the following equations: 

/ Afx JVx \ 


Ax Nt 


\j=0 
' N^ 


j=0 

Nt, 


Y1 Y.pZlk - 


( 2 ) 


i^l.Nf.k) +V'ij5 


n=0 m—0 


\k^0 


fe=0 


i = 0,...,Na;;j = 0,...,Nt, 

(5.38) 


where 


— '4’ (^x[Jj^ ^, tZNt.j) ’ * — 0 ,..., — 0 ,..., Nt- 


(5.39) 


Since Eormula (4.4) is exact for all polynomials hn{x) S P„, Equations (5.38) provide exact formulae for the bivariate 
shifted Gegenbauer interpolant, ^NtU, at the SGG nodes, , * = 0,..., iV^;; j = 0,..., N*. 


5.1. Global Collocation Matrix and Right Hand Side Constructions for Solving the Collocation Equations 

To put the pointwise representation of the linear system (5.36) into the standard matrix system form Ax = b, 
we introduce the mapping n = index(i, j) : n = i+ j {Nr^ + 1). Thus the matrix elements of the global collocation 
matrix A can be calculated by the following equations: 


-^index(z,j) 

,index(s,/c) ~ 

= {pZIs - 


-(2) 

Pl,Nt + l,s 

0 

ij^^PrZk 

+ p-^pZiZ) 

o' 

II 

o"' 

II 












(5.40a) 

-^index('i,j) 

,index(fc,j') ~ 

-- {pZlk - 


-'(2) 

. Pi.Ax + l, 

fc) 


i+p^pZlj 

+ , k — 0,... 

5 ^xi k ^ 

(5.40b) 

-^index(z,j^ 

),index(2,/c) ~ 

= (pZli - 


-'(2) 

,P/,Ax + l, 



+p^pZi,k 

) -pZi.k^ k = o, 

>...,Nt;k^j; 

(5.40c) 

-^index(i,j; 

),mdex(z,j') ~ 

■ {pZZ - 


-'(2) 

0 


+p^pZL 

+ i) -pZZj- 


(5.40d) 
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Algorithm A.l implements these formulas, and computes the right hand side of the linear system (5.36) by 
arranging the two-dimensional array (RHS)i,j = in the form of a vector array, {(RHS)„}^^q. Clearly, 

the construction of the global collocation matrix A requires the storage of its (L + 1) x (L -I- 1) elements, where 
L = Nx + Nt + NxNt.lt can be shown that the construction of the matrix A requires exactly 

1 + (1 + A,) (l + 5(1 + Ntf (1 + A,)) , 


multiplications and divisions, and 

2 + (1 + At) (1 + Aa;) (6 + 5Nx + At (5 + dA^,)), 
additions and subtractions, for a total number of 

A2 (A, + 1) (9A, + 10) + At (A, + 1) (19A, + 21) + 2A, (5A, + 11) + 15, 
flops. On the other hand, the construction of the right hand side requires 

(5 + 3Mt) (1 + At) (1 + Nx ), 

multiplications, and 3 additions and subtractions, for a total number of 

3 + (5 + 3Mt) (1 + At) (1 + Nx), 

flops. Hence, the total computational cost (TCC) of Algorithm A. 1 is given by 

TCC = 3Mt (At + 1) {Nx + 1) + A2 (A, + 1) (9A, + 10) + A* (A, + 1) (19A, + 26) + A, (lOA, + 27) + 23 
= O {NxNt{NxNt + Mt)), as Nx, Nt, Mt —>■ oo. 


For systems of small or moderate size, the solution by a direct solver is easy to implement; however, for large grids, the 
storage requirement of the global collocation matrix elements could be prohibitive, making fast iterative solvers more 
appropriate. Fortunately, the present numerical scheme converges exponentially fast for sufficiently smooth solutions 
using relatively small number of grids as we show later in Sections 5.3 and 6. In Section 6 we show also through a 
numerical example that the time complexity required for the calculation of the approximate solution at the 

collocation points ■, , i = 0,..., A^,; j = 0,..., Nt, using a direct solver implementing Algorithm A.l 

is approximately of O ((L + 1)^), as L —>^ oo, for relatively small values of Mt. 


5.2. Global Approximate Solution Over the Whole Solution Domain 

To approximate the unknown solution u at any point {x, t) s Df^, through Equation (5.37), we need to calculate the 
integrals Ji^x c\°'^{x),n = 0,..., Nx. Integrating Equations (A.11) in [23] on [—1,1], yield the following equations: 




Co°\cFo)d(To da I 



c[°\ao) dao dai 


= \{^ + xf = + C^^\x) + 4 “^ (c^2°'\x) - C^o‘^\x)'j , 

= ^ (-2 + x) (1 -f xf = (^-3 (a + 2) {2a + 1) C^o‘^\x) 

-3 (a + 1) (2 q; + 3) c[°'\x) + (a -f 1) (2a + 1) C';^“^(a;) + {a+ 2) {2a 


(5.41) 


1)), (5.42) 


^{-2 +x)x{l +xf, 
6 


j = 2 A a = 0, 


f-i III ('^o) dao dai = < 


1 


1 


4 (a + j) \a + j -I- 1 
j>2:j7^2Va7^0 


C'i+2(^) 


41 c‘°’(")) + <’ + 97’ (*)) , 
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where 


— (-9To(x) - 8(2ri(a;) + T 2 {x)) + T^ix )), j = 2 A a = 0, 
4o 


= < 


1 


1 


4(a+i) \a + j + l 
j >2 : j ^2y a^O, 




(5.43) 


c>(“) 


{x) 


,(a) 1 + 2 a 

=ATT^' 

j(a) _ ( 2 a + j )2 
“ (J+ 1)2 ’ 
,(a) _ 2 (g + j) 

a + j — 1 ’ 


j(“) _ _ (j 4)2 _ 

“A, (^ + _^-_i)(2„ + j-2)2’ 


j ^ 2 V a 7 ^ 0, 


4(-l)^ (2a-l) (g + j) 
(j +1)2 + j ~ 2)2 


((j^ + 2a (j + 2) - 4) a; + + 2a (j + 1) - l) , 


(5.44a) 

(5.44b) 

(5.44c) 

(5.44d) 

(5.44e) 


and {x)n = r(x + n)/r(x) = a: (x + 1) ... (x + n — 1), is the Pochhammer symbol (rising factorial). Therefore, the 
double integrals of cl^\x),j = 0,1,..., on [0, (], can be calculated as follows: 


(5.45) 


('^o) dao dcTi = ^ (^4“^ (x) + ( 2 ^) + ’ 

^ dao dai = ^x^ (-3^ + 2x) = (-3 (a + 2) (2a + 1) cl“\x) - 3 (a + 1) (2a + 3) (x) 


+{a + 1 ) ( 2 a + 1 ) 4 “^ (x) + (a + 2 ) ( 2 a - 1 )) , 


(5.46) 


Jo lo" dao dai = I 


— (Z - 2x) (3( - 2x) x^, j = 2 A a = 0, 


I 


^ { a + ] + l ^£+2^^) + C'g^(x)) + dg Cg4(x) + 9|J(x)^ , 


16 (a 

i>2:jV2Va7^0 


= < 


(-9Ti,o(x) - 8 (2r;_i(x) + Ti^ 2 {x)) + Ti 4 (x)), j = 2 A a = 0, 




16 (a + j) \a + j + 1 
t J > 2 : j 7 ^ 2 V a 7 ^ 0 


Jin; ‘^£+2(^)+4j + 45 ^.“^(a^)+^15 (a^)) > 


(5.47) 
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where Ti^„(x), is the nth-degree shifted Chebyshev polynomial of the first kind, and = 3^“^(2a;/Z — 1) Vj. 

Hence, 




^x{x-l) = ^4“^ 

^ 2x) [l-x)x = 

6 / 





P{l + 2a) 
48 (2 -h a) 




(5.48) 

(5.49) 


Ji,.c\f{x) = 


—- (I — x)x + Alx — 4x^) , j = 2 A a = 0, 

AA2(*) - (4j (45 + 45 <^4+2(^)) + > 

j>2:j^2VQ!^0 


= 


— {-9Tifi{x) - 8 Ti^2{x) +Ti^4{x) -I- 16), j = 2 A a = 0, 

45r (^44(2^) - 45 (45 (45 *^452;) + 45 ‘^£+2(^)) + ’ 

[j > 2 : jV 2 V a ^ 0, 


(5.50) 


where 


(a) 

, A- = 


(2 


lAa i6(a+j_l)2(2a + j-2)2^ 
1^2,j = {j - 1)2, 


(a) 

v\ ’ = 


(i + 1)2 (« + J +1)’ 

45 = (2a + j-2)2, 

45 = 2 (j + 1)2 (a -f j), 

45 = - (a + j - 1) (2a + 4)2; 

45 (^) = y (4 (a - 2)a -f 3) (a -f j - 1)3 ((-1)^ (Z - x) -f a;) . 

Using Equation (2.3), we can also show without stating the proof that 


Ji,^cl%\x) = 


{A {a-2) a+ 3) I 
4 (j + 1)2 (2a + j — 2)2 


(5.51a) 

(5.51b) 

(5.51c) 

(5.5 Id) 
(5.51e) 
(5.5 If) 
(5.5 Ig) 


hFi ( -j - 2,j -f 2a- 2;a - 1 - -f (-1)^ {x - 1) - ^ 

5 ,13 

2^“^2’2- 


(5.52) 

Using the above formulae, the SGPM directly approximates the solution at any point in the range of integration; 
on the other hand, finite-difference schemes, for instance, must require a further step of interpolation. 


5.3. Convergence and Error Analysis 

The following theorem gives the bounds on the discrete shifted Gegenbauer coefficients ^n,m Vn, m. 

Theorem 5.1 

Let u{x, t) S {Ffr) ’ the solution of Problem V. Suppose also that u is interpolated by the shifted Gegenbauer 
polynomials at the SGG nodes, (a::|5,, C = 0, ■ ■ ■ ^ Ffx;j = 0,..., Nt, on the rectangular domain Df.^. 

Then the discrete shifted Gegenbauer coefficients ^n,m, given by Equations (5.26) are bounded by the following 
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inequalities: 




4 (n + a) (m + a) r(n + 2a) T{m + 2a) 

r 2 ( 2 a + l)r(n + l)r(m + l) 

|(n + a) (m + a)| 


f +a- 1 


f+a-1 


4 \{n + a) {m + a)\ F + a) F + a) 


F2(q; + 1) F F ^'n? m? (n + 2a)^ (m + 2a) 


Uxx\\i^tx> Q; > 0, 

II II n m ^ I 1 

2’I 


n m , 1 


2 F (+ a) \{n + a) {m + a)\ 
a2 {m + 2a)T{^) F(a) 
2 F (+ a) I (n + a) (to + a) I 
[ a2 ^n(n + 2a)F(^^) F(a) 


a — 1 


f+a-l 


Ilo°(d2^)> 2 ^ Y ^ ~2 ^ ^ 

1 


li°°(-DY)’ 2 ^ ^ 2 ^ ^ 2 ^ ^ 

(5.53) 


Moreover, as n, to —oo, the discrete shifted Gegenbauer coefficients 4>n,m, are asymptotically bounded by: 


(nTO)^“ ||M;^a;|lioo(£,2 j, a > 0, 


< 




Z?2 (nTn) II^oo^£)2 ^ a <. 0 , 


where 


(a) Tt'^2^~*°‘ --4“-2 


Gi = 


F2 (a+i) F2(a+1)’ 


for some constants ^4^“^, > 1, dependent on a, but independent of n and to. 

Proof 

Since 

Nt 


(5.54a) 

(5.54b) 

(5.55a) 

(5.55b) 


/^(“) 

l^T^n, 7 r, 


£ Y1 ^s,k l,rCi% ti%,k) 


s=0 k—0 

l,T 

Nt 


\^n,mJ ^N^,Nt,s,k Vs,k ''At./c J > 


where 


s=0 fc=0 


\ (c>^) _ \ ('^) V/tt, 777, 

'^n^m v/t,//t, 


^Nt:,Nt,s,k ~ ^Nt,,s ^Nt.k^ S — 0, . . . , iVa;; fc — 0, . . . , iV^; 
Cn“i(a;, f) = (a:) C*") (f) Vn, to. 


Then, 


^n,m ^ 


7rF2 (a + i) 

F2 (a + 1) 


^(a) 


||'na;a;|| ^cxj/" £)2 , u — 0, . . . , Nx'l kfi — 0, . . . , N^, D - [ 

L°°(D^) v t<x/ 


(5.56) 

(5.57) 

(5.58) 


■1,1] X [-1,1]. 
(5.59) 
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Hence Inequalities (5.53) follow directly from Inequality (5.59), Equations (2.6), (5.56), (5.58), and Lemma 4.1 . Now, 
since 


r(a; -I- 1) < ■\/2 TT X ^ 




1 

1620 


Va; > 0; 


(5.60) 


cf. [27], then we can easily show that 



N, 



< 


^2§ 2a e n 2a 1 + q,) ^^n+2a+l ( 

n!r(a+i) r(a + l)(n + 2a) ^ ^ V 


1 


(n -I- 2a + 1) sinh 


1 


n -I- 2a -I-1 


-+a 


Va > 0. 


1620 (n-f 2a-f 1)'’ 
(5.61) 


Hence, 


^2 23—4a ^ — 4a —2 , 


< . + + . („ + 2a + l)”+2»+4 („ + 2o + l)”+=“+4 


i! ml (a -I- I) T‘^{a -|- 1) (n -I- 2a) {m 2a) 


1 


- r +1 

1620 (n-f 2a-f 1)^ J V1620 (m-f 2a + 1) 


1 


-f 1 ) ( (n -I- 2a -f 1) sinh 


1 


n -I- 2a -I- 1 


--l-a 


{m + 2a + 1 ) sinh 


1 


m -I- 2a -I- 1 


-+a 


J^xxW 


(5.62) 


d[°'‘{ nm)‘^°‘\\uxx \\(^02 )Va>0, asn,TO-)-oo, 


which proves the asymptotic inequality (5.54a). Similarly, we can prove the asymptotic inequality (5.54b) using 
Equation (4.29c). □ 


The following two corollaries underline the cases when the solution u of Problem 7^ is a polynomial or an analytic 
function. 


Corollary 5.1 

Let r s Z+, u{x,t) G P^, in x be the solution of Problem V. Suppose also that u is interpolated by the shifted 
Gegenbauer polynomials at the SGG nodes, ,i = 0,..., N^-jj = 0,..., Nt, on the rectangular 

domain Df^. Then there exists a positive constant C, such that the discrete shifted Gegenbauer coefficients ^n,m. 
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given by Equations (5.26) are bounded by the following inequalities: 

n m 1 

GZ+A-- <a<0, 

^eZ+A|^Z+A-^<a 
^^Z+A|eZ+A-^<a 
(5.63) 

where Ci = ACjP, is a positive constant independent of n and m. Moreover, as n,m —> oo, the discrete shifted 
Gegenbauer coefficients are asymptotically bounded by: 


< 


4 Cl (n + a) {m + a)T{n + 2a) r(m + 2a) 


r 2 ( 2 a + l)r(n+l)r(m + l) 




Of > 0, 


C; I (n + a) (m + a) I 


f+0-1 


f +0-1 


2 / \ 2 
4 C; I (n + o) (to + o) I r (+ o) E (+ o) 


r2(o + 1) r r y to^ (n + 2o)^ (m + 2a)^ 

2 G/ r (+ o) I (n + o) (to + o) I 


o 2 yTO(m + 2 o)r(^^^) r(a) 

2 G/ r (+ o) I (n + o) (to + o) I 


o2 (n + 2o)r(^i±i) r(o) 


f +a- 1 


+ 0—1 


iI“IIl~(d,%)> 

II«IIl~(d2^)> 


where 


lkllLoo(£,2^), 

D2,i inm)°‘ ||M||ioo(£,2 j, 


o > 0, 

< o < 0, 


(5.64a) 

(5.64b) 


D 


(o) 

1,1 


CiD\ 


(q) 


D 


(o) 

2,1 


Cid;\ 


(5.65a) 

(5.65b) 


and given by Equations (5.55). 


Proof 

The corollary follows from the inverse inequality of differentiation for any polynomial U{x) G P,.; cf. [28, Inequality 
(9.5.4)]. □ 


Corollary 5.2 

Eet u{x, f), be the solution of Problem V. Suppose also that U{x) = u {x, argmax^^^^^ \uxx{x, f)|), is an analytic 
function on [0, 1], and that there exist constants p > 1, and C{1, p), such that for every fc > 0, 


[/(fc) 


L“[0.i] 


<C{l,p) 


pk 


(5.66) 


If u is interpolated by the shifted Gegenbauer polynomials at the SGG nodes, “ 0: • • • j = 

0,..., Nt, on the rectangular domain Df.^, then the discrete shifted Gegenbauer coefficients (j)n,m, given by Equations 

This is a preprint copy that has been accepted for publication in Numerical Methods for Partial Differential Equations. 


























HIGH-ORDER NUMERICAL SOLUTION OE THE TELEGRAPH EQUATION 


23 


(5.26) are bounded by the following inequalities: 

8C{1, p) (n + a) (m + a) r(n -|- 2a) T(rn + 2a) 
r2(2a -f 1) r(n -L 1) T{m + 1) ^ 

2C{l,p) \ {n + a) (to -I- a) I 


(pa)2 


f +a- 1 


a > 0, 
+ a — 1 


8C{l,p) |(n + a)(TO + a)| r(^i^+a) r(^^+a) 

- r2(Q, + 1) r (nM) r w? {n + 2a)^ {m + 2a) 

AC{l,p)T + a) I(n-I-a) (to-I- a)I 


(pa)^ ^TO (to -I- 2a) r r(a) 

4C'(Z, p) r -I-a) I(n-I-a) (to-I- a) I 


f+a-l 


f+a-l 


n m 1 

eZ+A-- <a<0, 

2 , Y A-- <a<0, 


”eZ+A|^Z+A-^<a<0, 


^^Z+A|eZ+A-^<a<0. 


(5.67) 


Moreover, as n,m^ oo, the discrete shifted Gegenbauer coefficients (j)n,m, are asymptotically bounded by: 


^n,m 

< ‘ 






where 


Di,ppinm)°‘, a>0, 
D 2 lpinm)’^, -i<a<0, 


dTIp = ^c{i,p)d[^\ 
D2lp = \c{lp)Dt\ 


(5.68a) 

(5.68b) 

(5.69a) 

(5.69b) 


and d''P ; , are as given by Equations (5.55). 

Theorem 5.1 and its corollaries show that the coefficients of the bivariate shifted Gegenbauer expansions are 
bounded for — 1/2 < a < 0, as n,TO — > oo, but their magnitudes may asymptotically grow very large for a > 0, 
breaking the stability of the numerical scheme. Moreover, if we denote the asymptotic bound on the coefficients 
4’n,m, by to, a > 0, and by B~^n, to, a < 0, then Theorem 5.1 also manifests that the coefficients of the 

bivariate shifted Gegenbauer expansions decay faster for negative a-values than for non-negative a-values in the 
sense that 

B+ 

^— T —T——r(nTO)^ asn, TO—>'Oo, (5.70) 

where a+ and a_ are the Gegenbauer parameters of non-negative and negative values. This may suggest at first 
glance that numerical discretizations at the SGG nodes are preferable for negative a-values. However, we shall prove 
in the sequel that the asymptotic truncation error is minimized in the Chebyshev norm exactly at a = 0; that is, when 
applying the shifted Chebyshev basis polynomials. 

Remark 5.1 

Collocations at positive and large values of the Gegenbauer parameter a as n, to —oo are not recommended as can be 
inferred from Theorem 5.1 and its corollaries. However, the potential large magnitudes of the expansion coefficients 
are not the only reason for the instability of the numerical scheme in such cases. In fact, recently Elgindy and Smith- 
Miles [23] have also pointed out that the Gegenbauer quadrature ‘may become sensitive to round-off errors for positive 
and large values of the parameter a due to the narrowing effect of the Gegenbauer weight function w^°‘\xy, which 
drives the quadrature to become more extrapolatory; cf. [23, p. 90]. 

The following theorem highlights the total truncation error of the SGPM. 
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Theorem 5.2 

Consider the integral formulation of the hyperbolic telegraph PDE (5.13), and let 

where Q; rp, are as defined by Equations (5.14) and (5.15). Also let 

4^i,j — 4^i,j i^j y ^ — 0; ■ • ■ 5 ^xy j — O5 ■ ■ • ? ^ty 


> r = o,..., iV,; j = O,..., tv*, 


(5.71) 

(5.72) 

(5.73) 


be the integral hyperbolic telegraph PDE (5.13) at the interior SGG nodes iy'^^rlit G = 0; • ■ ■ y^x', j = 

0,..., Nt, and its discretization using the SGPM presented in Section 5, respectively. Then the total truncation error 

ETic/)^,J)= c|)^,jy r = 0,..., TV, ; J = 0,..., iV*, (5.74) 

is bounded by 

2-N^-Nt-Mt-4: 


K 


K 


A 


iVx + 1 


(iV, + l)!(At + l)!(Mt + l)! 


)Mt 


K 


Mt + l 


{Mt + 1)! (2^‘+2 Z^==+i [l + 


K 


(a) 
At + 1 


(TVt + l)!M;^ 


(x) 


a 




L°°[0,/] 


I _Aft + l M) 


a 


(a) 

l,Nt + l 


L“[0,t] 


)N^ + 1 


K 


(a) 
-^x + l 


(TV, + 1)!M;’ 


it) 


M 


+ +^TArtj + ( 2 , \Pl\ |/32|) 

(iV, + l)!(TV* + l)! 


(x,t) 


fA^x + l 


I 9A^a;+A£+2 ATt + l xitx) 


L°°[0,i] 

|/?2|+m)‘)) +2mW 1^1 i)) 


-^AT^ + l x^Nt + l 


c 


where 


= max 

(x,t)eDf^ 

= max 

(x,t)eDf,^ 

= max 

(x,t)eDf^ 

= max 

^ {x,t)&Df^ 

= max 

(x,t)eDf 




D 


(Aft + l) 


(j){x,t) 




f{x,t) 


D 


{N^+Nt+2) 
x,t 


^(Xyt) 


(a) 

(5.75) 

(5.76a) 

(5.76b) 

(5.76c) 

(5.76d) 

(5.76e) 


tj{N^ + 1) = rjiNt + 1) ^ . j^iN^+Nt+2) ^ ^^"+^*+2 

^ iV +1 ^ t - Tvr. I 1 1 Xy ^ ^ ^r i i tvt. 


dx 


dt 


TVt + 1 5 ^X,t 


dx 


Nx+i 




a+i 


Moreover, as TV,, TV*, Mt —> oo, the total truncation error is asymptotically bounded by 

ci(iV, A^tMt)“+5 

(c3 TV“+^ + C 4 4 ^^ TV^+2^ + C5 4 ^* iV“+^ iVf‘+2 ) + ce (^) 


Mt 




(5.77a) 


Cl (A, Nt MtY 


“+■5 A—Nic—Nt AT—N^—a — 2j\^—Nt — a — 2 ax 2 


N. 


■Mt 


C2 A 


(“) „Nt.^Nt 


C3A(“) 




+C6 4“' 


(I) 


Mt 




+ C5A4 4^*6^^^^^ a 4+®) 

-i < a < 0, 


(5.77b) 
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for some positive constants Cs,s = 1,..., 6; > 1, fc = 1,2,3, independent of N^;, Nt; Mt. 


Proof 

Consider 


the 


discrete 


approximations 


’ T.JVt.j ’ rNa:.i ’ t.Nj.j \ / ’ t.JVj.j ' ^ ’ T.JVt.j ' ^ 

defined by Equations (5.34) at the SGG nodes, l i) G = 0,..., j = 0,..., A^t, and denote the trunca¬ 
tion errors associated with them by J (^^1*^'^'/’) >^Li 1 ^i,j ) Ei,j 

Eij (^l 2 ^ f'^, respectively. Then 


\ET{f^,3)\ < \E^,j (j<^)| + \Pi\ \e,^3 + |/32| \e^, 3 [ii*^Jf)\ + \e,,3 (4*V)| 

Inequality (5.75) follows then from the error formulae (4.5) and (4.23) by observing that 
2_Ar.-2;N. + l + 




'’yJ \ 1 


S. 




{No, + 1)! 




||^(“) II 

II bN^ + l||^oo[o,i]’ 


\E. 


2—Nt — 2^Nt-\-l ^ 


<t> 1 

{Nt + 1)! 

l^(“) 



c. 


lia) 


'L^t-|-l||j;^oo[Q^T.]’ 


(5.78) 

(5.79a) 

(5.79b) 


o——3 A^t + l (g^) j.(g) / 

^ 


n^At) 1 ^ 

L-(“) 

{Nt + 1)! 

K^N) 
^Nt + l 


1 (A^^ + 1)! 

Nx-\-l 



3,,, (4*) J 0)1 


< 


-Nt-4 Nt + l (a) /.(a) )||r'(“) II 

(Nt + iy. 


2m1‘^ -P 


(5.79c) 

2-Ar,^Ar,+i^(x.i)|U(a) 


l,Nx-\-l I 


Il°^[c 


(U'v.)| < (0 




it) 


(Mt + iy. 


ll^(a) II 

II '^’^*+uIl~[o.t]’ 


Ei 


\3 {j2^ I 
E^,j {i^2 ^/)I 


< 


< 



y^T.Nt.j) 

^<'l 

{Mt + 1)! 

l^(“) 

|-^Mt + l 

V“) 'i 


{Mt + 1)! 

l^(“) 

|-^Mt + l 



c. 




(iV. + 1)! 

(5.79d) 

(5.79e) 


(5.79f) 


C 




(5.79g) 


The asymptotic inequalities (5.77a) and (5.77b) follow by taking the limits of Inequality (5.75) when Nx, Nt, Mt —> oo, and 
using Lemmas 4.1 and 4.2. □ 

The following corollary manifests that the asymptotic truncation error is minimized in the Chebyshev norm when 
applying the shifted Chebyshev basis polynomials. 

Corollary 5.3 

Let u{x, t) S {Efx) ’ the solution of Problem V. Suppose also that u is interpolated by the shifted Gegenbauer 
polynomials at the SGG nodes, ^a;|^ v'^^rNt ^ ’ * “ 0,..., j = 0,..., Nt, on the rectangular domain Df.,.. Then 
the shifted Chebyshev basis is optimal in the Chebyshev norm as Nx,Nt;Mt —> oo. 
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Proof 

If we denote the asymptotic truncation error bounds for positive and negative a-values by E^{(l)ij) and 
respectively, then we find by comparing the asymptotic formulae (5.77a) and (5.77b) that 

\E+{fi^^)\ > yi,j. (5.80) 


That is, the asymptotic truncation error bound is bigger for discretizations at the SGG nodes j, i = 

0,..., Nj;-, j = 0,., Nt, with positive a-values, and the gap grows wider for larger positive a-values. However, at 
a = 0, we find that the asymptotic truncation error bound, E^{(l)ij), is smaller in magnitude than \Ef{(l)ij)\, since 


\ ~ \^Ti(t>i,j) \ =Ci4 


-N.-Nt-Mt ]\f-N.-2 Mt § ^Nt+2 ^ 4“‘ 


+ - l) C6 4^=“ e^* -f C2 4“‘ e^‘ + ® \/N^Nt 

4“^ -^) + (^i“^ “ l) > OVi, j. (5.81) 

Hence the shifted Chebyshev basis is optimal in the Chebyshev norm as Nx,Nt\Mt —> oo. 


□ 


Notice that Corollary 5.3 is only valid for large numbers of expansion terms; however, the shifted Chebyshev 
basis is not necessarily optimal for small/medium range of expansion terms, where other members of the shifted 
Gegenbauer family of polynomials could exhibit faster convergence rates as recently shown by Elgindy and Smith- 
Miles [29] and Elgindy and Smith-Miles [24]. The following corollary highlights the convergence order of the SGPM. 

Corollary 5.4 

Eet u{x, f) G i^fr) ’ the solution of Problem V. Suppose also that u is interpolated by the shifted Gegenbauer 
polynomials at the SGG nodes, ^ ) * = Oj ■ ■ • j j = 0,..., iV^, on the rectangular domain Then 

the total truncation error of the SGPM is of 

/ / / \ ^ ^max \ 

/ / 1 \ 2 Nmax \ 

^ I \ ^ ^ ^5_g2b) 


as N^, Nt, Mt oo, where /max = max{/, t), TVmax = maxjiVj,, iVmin = minjiVa;, Nt, MJ. 

Proof 

The proof follows directly from Inequalities (5.77a) and (5.77b). □ 

Corollary 5.4 shows that the total truncation error of the SGPM decays faster than any finite power of 1/A^min, for 
— 1/2 < a < 0, exhibiting an “infinite order” or “exponential” convergence. 


6. NUMERICAL EXPERIMENTS 

In this section, we apply the proposed SGPM on four well-studied test examples with known exact solutions in the 
literature. Comparisons with other competitive numerical schemes are presented to assess the accuracy and efficiency 
of the SGPM. The numerical experiments were conducted on a personal laptop equipped with an Intel(R) Core(TM) 
i7-2670QM CPU with 2.20GHz speed running on a Windows 7 64-bit operating system. The resulting algebraic 
linear system of equations were solved using MATLAB’s “mldivide” Algorithm provided with MATLAB V. R2013b 
(8.2.0.701). The optimal S-matrix was constructed using Algorithm 2.2 in [23] with M^ax = 20. The change of 
variable (4.17) was applied using e = 3eps, where eps « 2.2204 x 10“^®, is the machine’s floating-point relative 
accuracy. The solutions a* of the one-dimensional optimization problems (4.16) were obtained using MATLAB’s 
line search algorithm “fminbnd” with the termination tolerance “TolX” set at eps/2. All of the test examples were 
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discretized at the shifted Chebyshev-Gauss nodes ^%,j) S 7 = 0,..., defined by 

(2j - l)7r 


^(0) AO) 

‘'t. 


3) \ 

,Nt,J J 


1 (^ ( (2* - l)7r 

- 1 — cos --- 

2 V V 2n 


1 — cos 


’ 2 V V 2n 

In all of the numerical experiments, we report the — , P—, norms of the absolute error matrix, 


= l,...,iV. (6.1) 


£ = 


„ fjo) ,(o) \ _ p / (0) (0) \ \ 


defined by 


I!^’ll2 

Halloo 


= max 

3 


N 

/ . ^ ^T,N,j J 

2=1 


= VAmax i£^£), 


= max 


N 

/ , I “ T,N,jj 


Pn,nu 


Pn,nu 


(0) 


(0) 


AO) 

^r,N,j 


AO) 

^T,N,j 


) 

) 


i=i 

where denotes the transpose of £ and Amax is the maximum eigenvalue of £'^ £. We shall also give the 

L°°-norm of the absolute error, \u{x, t) — Pn,nu{x, t)|, defined by 

Pn N = - PN,Nu{x,t)\\r^,r,2 ) = max \u{x,t) - PN^Nu{x,t)\, 

' ^ 0^' {x,t)eDf^ 

and provide the average elapsed CPU time in 100 runs (AECPUT) taken by the SGPM to calculate the approximate 
solutions. Moreover, in each test example, we plot the absolute error function, 

EN,Nix,t) = \u{x,t) - PN,Nu{x,t)\, {x,t) € Df .^, (6.2) 

and sketch the exact solution against its approximation over the entire domain Df^. 


Example 1 Consider Problem V with I = t = f3i = P 2 = ^,fix,t) =x‘^ +t— l,gi{x) = x‘^,g 2 {x) = l,hi{t) = 
t; h 2 {t) = 1 +1. The exact solution of the problem is u{x,t) = x"^ +1. The plots of the exact solution, its 
approximation, and the absolute error function on Df ^ using N = Mt = i, are shown in Figure 1 . The 
l^ — , norms of the absolute error matrix £, the L°°-norm of the absolute error, and the AECPUT are shown 

in Table I. The plots and the numerical results demonstrate the power of the SGPM, showing fast computations with 
errors of very small magnitudes using relatively small number of expansion terms. Moreover, Table I manifests the 
ability of the SGPM to achieve higher order approximations by increasing the number of expansion terms in the 
optimal S-quadrature (i.e. the value of Mt) while preserving the same size of the global collocation matrix A; hence 
the dimension of the linear system (5.36). A plot of the elapsed CPU time to compute the approximate solution, 

Pn,nu, at the collocation points j, j) >*>7 = 0,..., A^, versus the total number of unknowns (L -L 1) using 

various values of N and Mt is shown in Figure 2, where we observe an approximate time complexity ofO((L + l)^), 
as L —>^ 00 , for relatively small values of Mt. 


Example 2 Consider Problem V with l = T = l,pi = 10,/32=24:,gi{x)=x^{x—iy,g2{x) = 2x^{x — 
1)^,= h 2 (f) = 0;/(x, f) = 4e^* (x — 1)^ (l2a::^ — 24x^ — 2x^ + 14a: — 3 ). The exact solution of the 
problem is u{x, t) = x^(x — l)"‘e^*. The plots of the exact solution, its approximation, and the absolute error function 
on Dh using N — Mt = 8, are shown in Figure 3. The l^ — , P—, norms of the absolute error matrix £, the L°°- 
norm of the absolute error, and the AECPUT are shown in Table II. A comparison between Mohanty’s finite difference 
method [15], Ding et al.’s non-polynomial spline method [21], and the SGPM is also shown in Table III. The plots 
and the numerical comparisons show the rapid convergence rates and the memory minimizing feature of the SGPM. 
For instance. Ding et al.’s non-polynomial spline method [21] requires the solution of a linear system of equations of 
order 119 x 119, versus 9x9, for the SGPM to establish the same order of accuracy. 
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Figure 1 . The numerical simulation of the SGPM on Example 1 . The figure shows the plots of the exact solution u{x,t) on 
Dll (upper left), its approximation P4 4u(a;, t) (upper right), the absolute error function i?4 4(4;, t) (lower left), and its values 
at the final time, i?4 4(2:, 1 ) (lower right). The optimal S-matrix was constructed using Mt = 4 , and the plots were generated 
using 100 linearly spaced nodes in the x- and f-directions from 0 to 1. 


Example 1 

N,Mt 

4,4 

6,6 

6,9 

8,8 

8,12 

l|£|li 

1.665 X 10“i5 

1.810 X 10-i'‘ 

4.552 X 10-16 

3.408 X lO-ii 

1.554 X 10-11 

I|£|[2 

1.127 X 10“i5 

1.448 X 10-i'‘ 

2.947 X 10-16 

2.457 X lO-ii 

1.095 X 10-11 

II^IL 

1.665 X 10“i5 

2.087 X 10-i'‘ 

4.663 X 10-16 

4.025 X lO-ii 

1.649 X 10-11 

IT'00 

1.332 X 10“i® 

3.553 X 10-15 

2.887 X 10-16 

4.219 X 10-16 

4.663 X 10-16 

AECPUT 

0.012s 

0.028s 

0.029s 

0.279s 

0.290s 


Table I. The norms of the absolute error matrix S, the L°°-norm of the absolute error, and the AECPUT of 

Example 1 . 


Example 2 

N,Mt 

00 

00 

10,10 

12,12 

14,14 


5.418 X 10-1* 

1.883 X 10-11 

8.119 X lO-ii 

1.076 X 10-11 

I|fil2 

3.759 X lO-s 

1.542 X 10-11 

6.271 X lO-ii 

8.545 X 10-16 

II^IL 

5.089 X 10-® 

2.018 X 10-11 

8.952 X lO-ii 

1.158 X 10-11 

rpoo 

^N,N 

1.420 X 10-® 

4.222 X 10-1® 

1.331 X 10-11 

1.697 X 10-16 

AECPUT 

0.279s 

0.553s 

1.001s 

1.257s 


Table II. The —norms of the absolute error matrix S, the L°°-norm of the absolute error, and the AECPUT of 

Example 2 . 
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Time to Solve vs. Total Number of Unknowns 



Figure 2. Plot of the CPU time versus the total number of unknowns {L -\- 1) in log-log scale using N = 8(2)80 (i.e. 
80 < 1/ < 6560), and Mt = 4(2)16. The CPU time curves eventually approach a straight line with an average asymptotic slope 
of approximately 2, for increasing values of L, indicating an approximate time complexity of O ((L + 1)^), as L —>• oo, for 

relatively small values of Mt . 



Figure 3. The numerical simulation of the SGPM on Example 2. The figure shows the plots of the exact solution u{x,t) on 
Dll (upper left), its approximation Pg 8 M(a:, t) (upper right), the absolute error function Pg 3 ( 3 ;, t) (lower left), and its values 
at the final time, Pg 3 ( 3 :, 1) (lower right). The optimal S-matrix was constructed using Mt = 8 , and the plots were generated 
using 100 linearly spaced nodes in the x- and f-directions from 0 to 1 . 


Example 3 Consider Problem V with I = t = 1, ^i = 12, P 2 = 4 ,51 (a;) = sin{x),g 2 {x) = 0,hi{t) = 0,h2{t) = 
sin(l) cos{t); f{x,t) = 4 (cos(f) — 3 sin(f)) sin(a;). The exact solution of the problem is u(x, f) = sin(a;) cos(f).The 
plots of the exact solution, its approximation, and the absolute error function on Dfi using N = Mt = 4:, are shown 
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Example 2 

X 

Mohanty’s method [15] 

Difference scheme (11), 77 = 5 

Ding et al.’s method [21] 

— 12 j ^2 — g 

N = Mt = 8 

Present method 

N = Mt = 10 N = Mt = 12 

N = Mt = 14 

0.2 

4.956 X 10-^ 

1.922 X 10-*’ 

1.220 X 10-“ 

3.939 X 10-’“ 

1.266 X 10-” 

1.560 X 10-’“ 

0.4 

2.393 X 10-^ 

2.815 X 10-’“ 

2.740 X 10-’“ 

4.527 X 10-’“ 

4.819 X 10-’“ 

9.021 X 10-’“ 

0.6 

2.393 X 10-^ 

2.815 X 10-’“ 

2.740 X 10-’“ 

4.526 X 10-’“ 

4.774 X 10-’“ 

1.013 X 10-’“ 

0.8 

4.956 X 10-^ 

1.922 X 10-“ 

1.220 X 10-“ 

3.939 X 10-’“ 

1.252 X 10-” 

1.443 X 10-’“ 


Table III. A comparison of Example 2 between Mohanty’s finite difference method [15], Ding et al.’s non-polynomial spline 
method [21], and the SGPM. The table lists the absolute errors at a; = 0.2, 0.4, 0.6; 0.8, and f = 1. The results of Mohanty’s 
method [15] and Ding et al.’s method [21] are exactly as quoted from Ref. [21]. 


in Figure 4. The , P— , l°°—norms of the absolute error matrix £, the L°°-norm of the absolute error, and the 

AECPUT are shown in Table IV. A comparison between Dosti and Nazemi’s quartic B-spline collocation method 
[19], Mittal and Bhatia’s cubic B-spline collocation method [20], and the SGPM is also shown in Table V. Table IV 
shows one of the advantageous ingredients of the SGPM: “the ability to achieve higher-order approximations while 
preserving the same dimension of the linear system (5.36)”. Table V shows the power of the presented scheme, which 
constructs higher-order approximations using as small as 5 expansion terms in both spatial and temporal directions. 



Figure 4. The numerical simulation of the SGPM on Example 3. The figure shows the plots of the exact solution u{x,t) on 
^ 1.1 (upper left), its approximation P 4 4 u(a;, t) (upper right), the absolute error function i ?4 4 ( 4 ;, t) (lower left), and its values 
at the final time, i ?4 4 ( 0 ;, 1) (lower right). The optimal S-matrix was constructed using Mt = 4, and the plots were generated 
using 100 linearly spaced nodes in the x- and f-directions from 0 to 1 . 


Example 4 Consider Problem V with Z = r = l,/3i = 20,132 = 25, ( 71 ( 0 ;) = sinh(a:), < 72 ( 2 :) = —2sinh(a;), hi(f) = 
0, h 2 {t) = sinh(l); f(x, t) = —12e“^‘ sinh(a:). The exact solution of the problem is u{x, t) = sinh(a:). The 
plots of the exact solution, its approximation, and the absolute error function on Dh using N = Mt = fnre. shown 
in Figure 5. The P— , P— , l°°—norms of the absolute error matrix £, the L°°-norm of the absolute error, and the 
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Example 3 


N,Mt 

4,4 

4,5 

4,6 

6,6 

ll^lli 

2.384 X 10^'* 

1.036 X 10^^ 

9.391 X 10-® 

2.061 X 10-'^ 

I|f|l2 

1.762 X 10-'* 

7.438 X 10-® 

6.855 X 10-5 

1.533 X 10-^ 

ll^^lioc 

2.460 X 10-'* 

1.053 X 10-4 

9.731 X 10-5 

2.696 X 10-'^ 

rioo 

^N,N 

1.834 X 

8.060 X 10-® 

7.348 X 10-5 

1.160 X 10-'^ 

AECPUT 

0.012s 

0.012s 

0.013s 

0.027s 


Table IV. The noms of the absolute error matrix E, the Z/°°-norm of the absolute error, and the AECPUT of 

Example 3. 


Example 3 

t 

Dosti and Nazemi’s method [19] 

Mittal and Bhatia’s method [20] 


Present method 



h = 0.005, At = 0.001 

/i = 0.005, At = 0.001 

II 

II 

JV = 4, Mt = 5 

N = 4, Mt = 6 

N = Mt = 6 

0.2 

2.425 X 10-5 

6.827 X 10-5 

1.458 X 10-5 

7.449 X 10-5 

7.114 X 10-5 

4.188 X 10-5 

0.4 

7.932 X 10-5 

1.494 X 10-4 

4.905 X 10-5 

2.030 X 10-5 

1.922 X 10-5 

3.340 X 10-5 

0.6 

1.210 X 10-4 

2.241 X 10-4 

1.942 X 10-5 

1.127 X 10-5 

1.571 X 10-5 

9.030 X 10-5 

0.8 

1.488 X 10-4 

2.898 X 10-4 

8.293 X 10-5 

3.369 X 10-5 

3.154 X 10-5 

2.172 X 10-5 

1 

1.646 X 10-4 

3.439 X 10-4 

1.834 X 10-4 

8.060 X 10-5 

7.348 X 10-5 

1.160 X lO-’’ 


Table V. A comparison of Example 3 between Dosti and Nazemi’s quartic B-spline collocation method [19], Mittal and Bhatia’s 
cubic B-spline collocation method [20], and the SGPM. The table lists the maximum absolute errors at t = 0.2, 0.4,..., 1. The 
results of Dosti and Nazemi’s method [19] and Mittal and Bhatia’s method [20] are exactly as quoted from Ref. [20]. 


AECPUT are shown in Table VI. A comparison between Mohanty’s finite difference method [15] and the Pandit et al. 
combined Crank-Nicolson finite difference and Haar wavelets numerical scheme [6], and the SGPM is also shown in 
Table VII, which shows the root mean square error (RMS) for the three numerical schemes. Both tables show the fast 
execution times, the exponential convergence, and the cost economization features of the SGPM. 


Example 4 


N,Mt 

4,4 

4,5 

4,6 

6,6 

ll^lli 

4.027 X 10-3 

3.382 X 10-3 

2.807 X 10-3 

1.132 X 10-3 

11^112 

3.178 X 10-3 

2.688 X 10-3 

2.218 X 10-3 

8.012 X 10-5 

II^IL 

4.053 X 10-3 

3.578 X 10-3 

2.933 X 10-3 

1.251 X 10-3 

poo 

^N,N 

1.834 X 10-4 

2.571 X 10-3 

2.087 X 10-3 

4.855 X 10-5 

AECPUT 

0.012s 

0.013s 

0.014s 

0.029s 


Table VI. The P—, norms of the absolute error matrix £, the L°°-norm of the absolute error, and the AECPUT of 

Example 4. 


7. FUTURE WORK 

The present SGPM assumes sufficient global smoothness of the solution, and generally uses single grids for 
discretization on the spatial and temporal domains. An interesting direction for future work could involve a study 
of composite shifted Gegenbauer grids and adaptivity to improve the convergence behavior of the numerical solver 
when dealing with nonsmooth problems. On the other hand, the numerical experiments conducted in Section 6 
demonstrate the rapid convergence and stability of the SGPM; nonetheless, further stability analysis may be required 
to theoretically prove the stability of the SGPM on a wide variety of problems. 
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Figure 5. The numerical simulation of the SGPM on Example 4. The figure shows the plots of the exact solution u{x,t) on 
Dll (upper left), its approximation P 4 4 u(a;, t) (upper right), the absolute error function i ?4 4 ( 4 ;, t) (lower left), and its values 
at the final time, i ?4 4 ( 2 :, 1) (lower right). The optimal S-matrix was constructed using Mt = 4, and the plots were generated 
using 100 linearly spaced nodes in the x- and f-directions from 0 to 1 . 


Example 4 

Mohanty’s method [15] 

(2M)/(Af)/RMS 

The Pandit et al. method [6] 

(2M)/(At)/RMS 

Present method 

(V)/(Mt)/RMS 

(16)/(l/16)/0.3603 X 10-2 
(32)/(l/32)/0.6834 x 
(64)/(l/64)/0.2420 x 

(16)/(1/16)/4.4971 X 10-2 
(32)/(l/32)/5.8771 x IO -2 
(64)/(l/64)/6.3473 x IO -2 

(4)/(4)/6.382 X 10-^ 
(4)/(5)/5.375 X 10-^ 
(4)/(6)/4.436 X 10-^ 
(6)/(6)/1.184 X 10-® 


Table VII. A comparison of Example 4 between Mohanty’s finite difference method [15] and the Pandit et al. combined Crank- 
Nicolson finite difference and Haar wavelets numerical scheme [ 6 ], and the SGPM. The first two columns of the table lists the 
results in the form (2M)/(A t) /RMS, and are exactly as quoted from Ref. [ 6 ]. The last column of the table lists the RMS errors 

of the SGPM associated with various values of N and Mt . 


8 . CONCLUSION 

In this work, we developed a novel SGPM for the solution of the telegraph equation provided with some initial 
and boundary conditions. The method recasts the problem into its integral formulation to take advantage of the 
stability and well-conditioning of numerical integral operators. The discretization is carried out using some novel 
shifted Gegenbauer integration matrices and optimal shifted Gegenbauer integration matrices in the sense of solving 
the one-dimensional optimization problems (4.16). We established Algorithm A.l for the efficient construction 
of the global collocation matrix and the right hand side of the resulting linear system, which together with a 
standard direct solver can produce very accurate approximations. The TCC of the developed algorithm scales like 
O {Nj;Nt{NocNt + Mt )), as Nj;, Nt, Mt —> 00 . A numerical study on the time complexity required for the calculation 
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of the approximate solution at the collocation points using a direct solver implementing Algorithm A.l shows 
that it scales like O [{L + 1)^), as L ^ oo, for relatively small values of Mt, where (L + 1) is the total number 
of unknowns in the resulting linear system. Theorem 5.1 and its corollaries demonstrate that the coefficients of 
the bivariate shifted Gegenbauer expansions decay faster for negative a-values than for non-negative a-values. 
In fact, we proved that the coefficients of the bivariate shifted Gegenbauer expansions are bounded for a non¬ 
positive Gegenbauer parameter, a < 0, as n, m —oo. Corollary 5.3 shows also that the asymptotic truncation error 
is minimized in the Chebyshev norm exactly when applying the shifted Chebyshev basis polynomials. Corollary 
5.4 proves the exponential convergence exhibited by the SGPM. The extensive numerical results and comparisons 
demonstrate the fast execution, the exponential convergence, and the computational cost effectiveness of the proposed 
method. Moreover, the results show the ability of the numerical scheme to achieve higher-order approximations while 
preserving the same number of the solution expansion terms; thus the dimension of the resulting linear system of 
algebraic equations (5.36). The method is memory minimizing, easily programmed, and can be efficiently applied 
and extended for the solution of various problems in many areas of science. 
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A. A COMPUTATIONAL ALGORITHM FOR THE CONSTRUCTIONS OF THE GEOBAE COEEOCATION 
MATRIX AND THE RIGHT HAND SIDE OF THE RESUETING COEEOCATION EQUATIONS 


Algorithm A.l Matrix and Right Hand Side Constructions for Solving the Collocation Equations (5.36) 

Input: Positive integers N^, Nt, and Mu Positive real numbers I and r; constants /3i and ^ 2 ; SGG nodes x\°‘^ i 
0,..., Nx', the S-matrices Pp\ Pr^\ and Pr^^; the optimal S-matrices P^^^ and P^^^. 

L •<— Nx + Nt + NxNt 

for m = 0 to L do 
for n = 0 to L do 

-^n,m ^ 0 

end for 
end for 

for i = 0 to Nx do 




„(“) 


for j = 0 to Nt do 

n ^ index (i,j) 

An n ^ (p\,% - 


(„) 


for fc = 0 to Nx do 
if fc ^ j then 

■^ri,index(fcj) 

end if 
end for 

for fc = 0 to A^t do 
if fc ^ j then 

.^n,index(2,fc) 

for s = 0 to Nx do 
s ^ i then 


-'( 2 ) 

. Pi,k+i 


,.) (A: 


i(l) 


■ 1^2 P 


( 2 ) 


( 2 ) 


M) 


-(2) 

. Pl,Nx + l 


,k) (/3ipSj 


-P- 




+ Prik + P2P^r,ik) - P?j,k 


^nandex(s,/L‘) 




Pys - 


4 A2) 

end if 
end for 
end if 
end for 

(RHS)„ ^ - (pi pS.feV'i.a.fc) + ^2 E Prik^^.u,k) 

\ k=0 k=0 

end for 
end for 

return A; {(RHS)„}Eo 


,s) + (^2Pxlk) 


^ ( 2 ) f 

2^ PT,j,k Ji,U,k) 
fc =0 
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